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Abstract 

Using an updated simulation tool, we examine molecular junctions comprised of benzene- 
1 ,4-dithiolate bonded between gold nanotips, focusing on the importance of environmental 
factors and inter-electrode distance on the formation and structure of bridged molecules. We 
investigate the complex relationship between monolayer density and tip separation, finding 
that the formation of multi-molecule junctions is favored at low monolayer density, while 
single-molecule junctions are favored at high density. We demonstrate that tip geometry and 
monolayer interactions, two factors that are often neglected in simulation, affect the bonding 
geometry and tilt angle of bridged molecules. We further show that the structures of bridged 
molecules at 298 and 77 K are similar. 

Keywords: Molecular Junction, Molecular Wire, Molecular Electronics, Molecular Simula- 
tion, Mechanically Controllable Break Junction, Electron Transport, Single-Molecule Conduc- 
tance, Gold Nanowire, Benzenedithiol. 

Conductance measurements through molecular junctions have been at the forefront of nanoscale 
research for over a decade.l^^This work is motivated by the potential for fabrication of molecular- 
based electronic circuit elements'^ and, perhaps more so, discrepancies in the experimentally!^ and 
theoreticallyl^ reported conductance through a single molecule. The discrepancies have improved 
over the years, ^ due in part to the development of highly automated and optimized experimen- 
tal techniques (e.g., scanning tunneling microscopy break junction method,!^^ nanofabricated 
mechanically-controllable break junction technique^^, as well as the emergence of theoretical 
tools (e.g., self-consistent GW calculations,^^ approximate self-interaction corrections^^^l^ ca- 
pable of more accurately describing the HOMO-LUMO gap and energy level lineup between a 
molecule and two leads. Moreover, it has been repeatedly demonstrated that a spectrum of struc- 
tures exist in the experiments, some of which seem to appear more frequently than others based 
on the relative peak heights in histograms of the conductance.'^^^^^^^^'^ For example, recent low- 
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temperature (4.2 K) measurements of benzene- 1,4-dithiolate (BDT) showed several peaks between 
10^-^ Go and 0.5Go, where Go=^.^^ Results such as these have shifted focus away from reproduc- 
ing a single value of conductance towards, more generally, determining the structures responsible 
for the most-probable conductance values in a given experimental setup. Taking cues from exper- 
iments, researchers on the theoretical side have recently begun calculating the conductance of an 
ensemble of molecular junction structures.'^^'^^'^ Structures are obtained using molecular dynam- 
ics simulations in which the molecular junction is evolved through mechanical elongation! ^ ^ 1 ^^^ 
or compression,!^ or by thermal activation. 1^^^ Valuable information about how local structural 
conformations {e.g, oligomeric gold-thiolate units^^ and tilt angle^^) influence the trends in con- 
ductance has been provided by these studies. However, environmental factors (e.g., monolayer 
interactions, non-ideal electrode geometry) have not yet been included in these simulations, de- 
spite the fact that they are likely to influence the results.!^^^ 

Balancing accuracy and computational efficiency is a major challenge for simulations of molec- 
ular junctions. Simulations need to accurately capture the preferred bonding geometries while also 
incorporating environmental factors found in experiment. Quantum mechanical (QM)-based meth- 
ods, such as density functional theory (DFT), are capable of accurately resolving molecular- level 
bonding, but the high computational cost of QM methods may limit the system size, reduce the 
total number of independent statepoints, and require simplifications to the local junction environ- 
ment (e.g., neglecting monolayer effects, employing ideal electrode geometries, and considering 
single-molecule junctions only) _[1512il t26 | 32H 3l Additionally, energy minimizations often included in 
DFT calculations^^^^^ may produce configurations that are not likely for thermal systems. Meth- 
ods based on classical force fields have also been used to simulate molecular junctions'^^'^ and 
related systems.l^^^^ Classical force field (CFF) methods (i.e., molecular dynamics - MD - and 
Monte Carlo - MC - simulation ■^^) are able to handle larger system sizes and more statepoints 
than QM methods; however, metal- molecule interfaces exhibit complex bonding with preferred 
bonding sites that cannot be easily captured by conventional CFF models and methods. Previous 
CFF-based MD simulations of molecular junctions have only considered ideal junction environ- 
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ments, e.g., a single molecule sandwiched between perfectly flat electrode surfaces.'^^'^^'^ This 
is in contrast to experimental systems, where the bridged molecule may be surrounded by other 
adsorbed molecules (i.e., a monolayer) with electrodes that have curved geometries resulting from, 
e.g., the rupturing of a Au nanowire (NW), as carried out in the mechanically-controllable break 
junction (MCBJ) experimental technique .'^'^^^ Reactive force fields (e.g., ReaxFF) have shown 
promise as a compromise between QM and CFF methods,!^ but parameters for metal-molecule 
systems are still under development. 

Here, we present an updated CFF method that balances computational efficiency and accuracy. 
The method builds from our previously developed hybrid MD-MC methodic and is capable of 
incorporating conditions more representative of experiment than previous work. The technique 
allows larger systems with multiple bridged molecules, a large number of statepoints, and more 
realistic environmental factors to be considered, while retaining high accuracy through the use of 
bonding parameters derived from DFT calculations.'^ Additionally, we extend our prior workl^ 
by performing simulations within the semigrand canonical ensemble,!^ which includes MC moves 
designed to improve sampling of the preferred metal-molecule bonding geometries, further differ- 
entiating our technique from previous CFF-based studies. Using this approach, we are able to sim- 
ulate important aspects of MCBJ experiments (see Figure 1), from formation of a self-assembled 
monolayer (SAM) onto a Au NW surface, to elongation and rupture of the NW, and finally to 
trapping of a small number of molecules (between 1 and 4) in a break junction. 

We demonstrate the novelty and utility of the new simulation technique by generating a sta- 
tistical ensemble of Au-BDT-Au junctions, examining in detail the number, tilt angle, and bond- 
ing geometry of bridged molecules. By varying the extent of monolayer coverage, we find that 
monolayer packing is influential in the formation of single and multi-molecule junctions. We also 
perform simulations with idealized tip geometries both with and without monolayers, to allow 
comparison with systems commonly used in simulation and theory.^^ ^^ l "^^* ^ Lastly, we demon- 
strate that temperature can be used to control the number of bridged molecules. The computational 
tractability of the simulation method allows us to perform over 1,000 simulations, resulting in 
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Figure 1: Simulation snapshots of the MCBJ method, (a) BDT self-assembles onto an unstretched 
Au NW; a closeup is shown in (b). (c) Au point contact in the necked region of the NW after ~3.5 
nm of elongation, (d) Following NW rupture, the bulk BDT is evaporated from the simulation box. 
(e) The ruptured NW tips are brought together, resulting in the formation of a molecular junction. 

statistics on par with experiment. 

Results 

To generate non-ideal electrodes that are representative of those found in MCBJ experiments, we 
first perform ten independent simulations of the elongation and rupture of BDT-coated 1 .9-nm- 
diameter Au NWs (see Methods). The NWs are elongated in the [001] direction at a rate of 1 m/s 
and temperature of 298 K using a hybrid MD-MC technique. The next step in the MCBJ process, 
and the aspect we focus on in this paper, is the formation of a molecular junction, which we 
simulate using a MC-based method. Coupling each ruptured NW tip with one another (including 
a tip with itself) yields a total of 210 unique electrode-electrode combinations for performing 
simulations of the molecular junction formation process. We simulate a mixture of two BDT 
species, one of which bonds at on-top sites while the other bonds at on-bridge sites. Previous 
experimental'^^'^ and theoretica P^ I ^^^ ^ studies have demonstrated that the on-bridge site is the 
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energetically favored bonding site for benzenethiolatel^SMIlZl and alkanethiolates,'^^'^^'^ while on- 
top sites are important in low-coordination environments.'^ We implement BDT identity swap 
moves during MC sampling, which allows BDT to chemisorb to preferred sites. The interaction 
potentials for BDT bonded at on-top and on-bridge sites were previously fit to DFT calculationSj'^Sl 
thus enabling us to perform simulations in a classical framework while retaining high accuracy. 

Following NW rupture, we displace the Au tips in the x-y plane such that the bottom-most and 
top-most Au atoms in the top and bottom tips, respectively, are aligned along the z axis. Next, 
the tips are gradually pushed together, from Z = 20AtoZ = 6A (where Z is the inter-electrode 
distance), over the course of 25 million MC moves. We end each run at Z = 6 A since direct 
tunneling between electrodes has been shown to occur for Z < 6 A.^^^l^ Figure 2a shows a typical 
plot of the number of bridged BDT molecules as Z is decreased. Initially, at large values of Z, zero 
molecules are chemically attached to both electrodes. At Z < 1 1 A, a single-molecule junction 
forms, as shown in Figure 2b. As Z is decreased further, two (Figure 2c) and eventually three 
(Figure 2d) molecules connect in parallel. All images in this paper were rendered using Visual 
Molecular Dynamics.'^ 

Surface Coverage Effects 

We first explore the impact of monolayer packing by performing 210 simulations for each of four 
different surface coverages: 0.30, 0.40, 0.50, and 0.65 +1- 0.03. Surface coverage is defined here 
as the number of adsorbed molecules divided by the number of Au surface atoms. 0.65 +1- 0.03 is 
the maximum surface coverage obtained for the 20 ruptured Au NW tips, which closely matches 
the reported coverage for alkanethiolates on Au nanoparticles of diameter 1.3-1.4 nm.!^ 

Using molecule number data such as those shown in Figure 2a, we construct histograms (see 
Figure 2e) of the number of bridged molecules as a function of Z, with separate panels representing 
(from top to bottom) decreasing surface coverage and the color of the histogram bars corresponding 
to the number of bridged BDT molecules {n). We observe that the histograms of bridged molecules 
tend to increase with decreasing Z, with the exception of the n = I case, which exhibits a peak at 
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Figure 2: (a) Typical plot showing the number of bridged BDTs as the inter-electrode separation, 
Z, is decreased. This particular simulation results in (b) one bridged molecule from Z ~ 10.6-8.4 
A, (c) two bridged molecules from Z ~ 8.4-7.6 A, and (d) three bridged molecules from Z ~ 
7.6-6.0 A, with the corresponding images shown below. The bridged and non-bridged BDT are 
rendered differently in the images for clarity, (e) Histograms of the number of bridged molecules 
as a function of Z. The histogram bar colors correspond to the number of bridged molecules. The 
red arrows indicate the maximum Z at which the single-molecule histograms are at least 98% of 
their peak values. 

aU four surface coverages. These peaks, which are indicated with red arrows, appear due to the 
rate of formation of multi-molecule junctions exceeding that of single-molecule junctions; these 
peaks shift to higher Z for lower surface coverages. 

We also observe in Figure 2e that n depends on surface coverage. For most Z, the formation of 
at least one bridged molecule (n > 0) is most likely for surface coverage 0.50 and least likely for 
0.30. The optimal surface coverage for forming a single bridged BDT (n = 1) depends on Z; for 
Z > 10 A, intermediate coverages (0.40/0.50) provide the highest probability, while for Z < 9 A, 
n = 1 is most probable at maximum coverage (0.65 +/- 0.03). Low surface coverages (0.30/0.40) 
tend to result in the highest occurrence of multi-molecule junctions. Experimentally, conductance 
histograms often exhibit peaks at integer multiples (n) of a fundamental conductance value, with n 
corresponding to the number of molecules in the junction.l ^ l ^^l^^ ^^Two- and three-molecule peaks 
often occur in break junction experiments ;'^^^2IH four-molecule peaks have also been observed.^^ 
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These data match our results. Additionally, the relative peak heights in experiment generally de- 
crease with larger n, which from Figure 2e holds for most surface coverages and values of Z in 
our simulations. We conclude that the trends we observe in our simulations are in good agreement 
with experimental results, thus validating our methodology. 

It is important to note that surface coverage generally varies between experimental setups, 
with some experiments conducted at low coverages in order to provide available bonding sites 
for molecular bridging J 1 1 1 20 | 2 1 1 5 1 152| others performed with the bridging molecules diluted in 
a dense matrix of non-bridging adsorbate molecules .1^^^ In the seminal work of Reed and co- 
workers,'^the break junction was exposed to a solution of BDT for a long period of time, resulting 
in a densely packed monolayer on each of the Au nanotips. Subsequent theoretical work^S! sug- 
gested that the low conductance observed by Reed and co-workers could be attributed to weak 
electrical coupling between two overlapping BDT molecules; in other words, chemical contact be- 
tween a single molecule and the two electrodes was not established, owing to the lack of available 
bonding sites on each nanotip. Our results show evidence of such effects, but not to the degree that 
a single-molecule junction cannot form. That is, we observe that squeezing a single molecule into 
an already dense monolayer is compensated by the addition of a S-Au chemical bond; however, 
the energetic penalty for fitting more than one molecule is often too great to overcome. Note that 
the tip curvatures in our simulations may differ from those of Reed et al.,^ which may influence 
whether a molecule is able to bridge in densely packed monolayers. We do not consider such 
effects here. 

In addition to changing the number of available bonding sites, the packing density of a mono- 
layer also affects the mobility of adsorbed BDT, and thus influences whether a molecule can adopt 
one of the specific geometries required for bridging. The reduced interactions between adsorbed 
BDTs, along with the increased availability of bonding sites, is the cause of the shifting in single- 
molecule peaks to larger Z at lower coverages, as a second molecule can more easily bridge. This 
is also the cause for the large n>0 histograms at intermediate coverages, and large multi-molecule 
histograms at low coverage. It is somewhat surprising that the formation of three or four bridged 
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BDTs is more likely at low coverage since one might expect the number of molecules on each 
tip to be the dominant factor in determining the number of bridged molecules.'^ We point to the 
reduced monolayer interactions as the cause for this somewhat counterintuitive behavior. We also 
note that in experiments conducted at low coverages, there is often evidence of multi-molecule 
junctions J i ^l^o i ^ 1 1 While the exact surface coverage in these experiments is unknown, Figure 2e in- 
dicates that the relative frequency at which multi-molecule junctions form will depend on Z and 
surface coverage. 

Role of Non-Ideality 

In order to examine the impact of realistic environmental features, we next compare our MCBJ 
simulation results with those for idealized systems. First we explore the effect of using an ideal tip 
geometry. The ideal tip, which is shown in Figure 3a, is an atomically sharp, pyramidal structure, 
reminiscent of the electrode geometries used in numerous previous theoretical studies .[I Bl22|42t g4l 
Note that for the remainder of this paper we will only consider intermediate surface coverages of 
0.40. 

Figure 3b plots histograms of the number of bridged BDT molecules as a function of Z, with 
the ideal and NW tip results shown at top and bottom, respectively. The histograms demonstrate a 
tip geometry dependence; the probability of having n > is higher for the ideal tips at Z < 10 A, 
while the ideal tip histograms change more rapidly than those for the ruptured NW tips. We further 
assess the impact of tip geometry by plotting the bonding geometry as a function of Z, shown in 
Figure 3c. The separate panels display the three possible combinations of sites (i.e., on-top/on- 
top, on-top/on-bridge, and on-bridge/on-bridge) binding a bridged molecule. In general, on-bridge 
sites become more available for molecular bridging at lower values of Z, especially for the ideal 
tips where only on-top sites are accessible for bridging at high Z. In contrast to the ideal tip, a 
ruptured NW tip can be relatively flat at its apex, with on-bridge sites accessible at high Z. Lastly, 
we plot the tilt angle of bridged molecules, shown in Figure 3d. Before discussing these results, 
we introduce a simple compression model as a first approximation for relating the inter-electrode 
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Figure 3: (a) Molecular junction composed of two ideal, atomically sharp Au tips (at Z = 10.15 
A) and a single bridged BDT molecule. The bridged BDT and non-bridged BDTs are rendered 
differently for clarity, (b) Histograms of the number of bridged molecules at various values of Z. 
(c) The bonding geometry for bridged BDT molecules as a function of Z. Each panel represents 
the fraction of different combinations of on-top and on-bridge bonding, with (from top to bottom) 
"T-T" denoting on-top bonding at both tips, "T-B" denoting on-top bonding at one tip and on- 
bridge bonding at the other, and "B-B" denoting on-bridge bonding at both tips, (d) The tilt angle, 
6, versus Z. 

separation, Z, to the tilt angle, 0: 

Z(0) = Ds-scosO +2^Dl_^^-Dl_ssin^d, (1) 

o 

where Ds s is the distance between S atoms in a BDT molecule (6.28 A for our rigid model of 
BDT) and Ds au is the equilibrium S-Au bond distance (2.29 A for on-top bonding). This model 
assumes that the S atoms remain bonded to the on-top sites of each tip (with Ds au = 2.29 A), the 
BDT center-of-mass falls along the z axis made by the two Au tips, and the tips are aligned in the 
x-y plane. We find that these first two assumptions often break down for low Z; nonetheless, eq. 
1 establishes a baseline for comparison of tilt angle data, and qualitatively captures the behavior 
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expected from a bridged molecule that remains at the tip apex while compressed, as opposed to 
one that migrates to sites along the side of a tip. We note reasonable agreement between our 
compression model and the tilt angle data in Figure 3d, especially for the ideal tips. While the tilt 
angle trajectory of any single bridged molecule may differ significantly from the average behavior, 
as evidenced by the large uncertainty bars, the average trends are in qualitative agreement with the 
compression model, suggesting that molecules tend not to migrate to sites along the sides of the 
tips. For Z < 10 A the non-ideal tips result in tilt angles that are, on average, higher than those for 
ideal tips, indicating that the migration of bridged molecules to sites along the side of the tips is 
less common in systems with non-ideal tips. 

We next demonstrate the effect of a monolayer on the bonding geometry and tilt angle of 
bridged molecules. After obtaining twenty different monolayer arrangements on the ideal tip (Fig- 
ure 3a), we perform 210 simulations using each unique combination of the twenty BDT-decorated 
tips. We identify twelve runs resulting in the formation of a single-molecule junction at Z > 11 

o 

A. Using these single-molecule structures as the starting point, we then perform simulations in 
which the remaining monolayer molecules are absent from the electrodes, enabling us to assess 
the impact of adsorbate interactions with the bridged molecule. 

In Figure 4a we present the bonding geometry of the bridged molecules, with similar trends 
observed for the monolayer and no-monolayer scenarios, but quantitative differences. Recall that 
high monolayer density limits the availability of bonding sites while also reducing molecular mo- 
bility, which is responsible for the larger on-bridge peak in the no-monolayer systems shown in 
Figure 4a. To further investigate why the bonding geometry changes with Z, it is instructive to ana- 
lyze the S-Au bond energy. In Figure 4b we present the average S-Au bond energy versus Z for the 
no-monolayer runs. Because there is no monolayer present, the bridged molecule is able to freely 
explore the energetically favored sites at each tip. At large Z, molecular bridging is only possible 
with on-top/on-top bonding geometry; as Z is decreased, the energetically more stable on-bridge 
sites become accessible for bridging; for low values of Z, the compression of the tips gives rise 
to situations where on-bridge/on-top bonding geometry becomes energetically competitive with a 
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somewhat strained on-bridge/on-bridge connection. 
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Figure 4: Comparison of molecular break junction behavior in presence and absence of a mono- 
layer, (a) The bonding geometry for bridged BDT molecules as a function of Z. See the caption 
in Figure 3 for definitions of the abbreviated terms, (b) S-Au bond energy versus Z for a single 
bridged BDT molecule, neglecting monolayer effects, (c) The tilt angle, 0, as a function of Z. The 
compression model (eq. 1) is plotted as the red curve. The inset histograms show the distribution 
of tilt angles. 



Figure 4c plots the tilt angle versus Z of bridged molecules in the presence and absence of a 
monolayer. The compression model (eq. 1) is also shown (red curve) for reference. We find that 
the tilt angles of bridged BDTs for Z > 9.5 A agree closely for the cases where a monolayer is 
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present and absent. This regime is characterized by low tilt angles and, for Z > 1 1 A, long S-Au 
bond lengths. The maximum value of Z for which a bridged molecule forms is 12.11 A. This 
value of Z requires an average S-Au bond length of 2.92 A, in close agreement with the reported 
S-Au bond rupture distance of 2.86 A.I^For Z < 9.5 A, the monolayer and no-monolayer results 
differ markedly. In the presence of a monolayer the tilt angles of bridged molecules trend upward, 
indicative of the confinement of bridged molecules to the tip apex. In absence of a monolayer, the 
bridged molecules exhibit different tilt angle behavior, undergoing abrupt changes that coincide 
with changes in the bonding geometry (see Figure 4a). The inset in Figure 4c shows the entire 
distribution of tilt angles of bridged molecules. Bridged molecules reach a maximum of ~30° in 
absence of a monolayer, and exhibit two preferred tilt angles at 2.5° and 14.5°. On the other hand, 
the tilt angle distribution for bridged molecules in the presence of a monolayer is relatively flat 
from 0=2-35°, with a maximum value of ~50°. 

The differences we highlight for idealized systems are significant since the bonding geometry 
and tilt angle of bridged molecules have been demonstrated to affect experimentally observed prop- 
erties, namely conductance and inelastic electron tunneling spectra (lETS). Conductance has been 
shown to scale linearly with the number of bridged molecules,!^ while various studies^ ^'^ l ^^ l ^^ l 
have demonstrated that bonding geometry and tilt angle can affect conductance by an order of 
magnitude or more. For example, Haiss and co-workers showed that increasing the BDT tilt 
angle from = 0° to = 50° results in close to an order of magnitude increase in conductance, 
with the most pronounced increases occurring between 6 = 30-50°. Recall from the histograms 
in Figure 4c that the maximum tilt angle with a monolayer present is ~50°, but only ~30° with 
no monolayer. Thus, in this case, neglecting monolayer effects could result in significant under- 
predictions of conductance. In addition to affecting conductance, bonding geometry and tilt angle 
have also been shown to influence the lETS of molecular junctions.!^ 
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Figure 5: Comparison of molecular break junction behavior at 298 and 77 K. (a) Histograms of the 
number of bridged molecules at various values of Z. (b) The bonding geometry for bridged BDT 
molecules versus Z. See the caption in Figure 3 for definitions of the abbreviated terms, (c) The 
tilt angle, 0, as a function of Z. 

Role of Temperature 



The results presented until now have been for a temperature of 298 K. We now consider a tem- 
perature of 77 K, which corresponds to cryogenic conditions and has been used in experiments 
of Au-BDT-Au junctions.'^ We employ the same twenty ruptured NW tips for both temperatures, 
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performing 210 simulations in each case. Figure 5a shows histograms of the number of bridged 
molecules as a function of Z, at 298 and 77 K. Clearly, 298 K results in a significantly higher 
probability of forming a molecular junction composed of any number of molecules, for a majority 
of Z; thus, the reduced mobility of the BDT molecules at 77 K is detrimental to molecular bridg- 
ing. We further examine the influence of temperature by plotting the bonding geometry and tilt 
angle in Figure 5b and Figure 5c, respectively. Overall, the quantitative differences between the 
two temperatures are small. In Figure 5b, the fraction of on-top sites binding a bridged molecule 
is slightly higher at 77 K and low Z. We attribute this to molecules being unable to migrate off 
of on-top sites after bridging there at large Z. This explanation is corroborated by data in Figure 
5c, which displays higher tilt angles at 77 K than 298 K, indicating that confinement of bridged 
molecules to the tip apex takes place more often at the lower temperature. 

To our knowledge, no experimental data has been reported comparing the number, bonding 
geometry, or tilt angle of bridged molecules at different temperatures. Studies on temperature- 
dependent behavior have focused on other properties such as mechanical stabihty^^ and conduc- 
tance.l^^MCBJ studies of Au-BDT-Au junctions at 77 I^'^and 4.2 Kl^^have shown discernible 
peaks in histograms of conductance, but in neither case was an analysis of the relative peak heights 
at different temperatures reported. At 77 K, Tsutsui and co-workers'^ observed a peak in the BDT 
conductance histogram at 0.01 1 Go, matching the reported value at 298 K;'^^^this finding not only 
implies that coherent tunneling remains the dominant electron transport mechanism in the temper- 
ature range, but also that the most frequently occurring structures at 77 K and 298 K are similar. 
Our results in Figure 5b and Figure 5c support this latter conclusion, especially for high Z, as the 
bonding geometry and tilt angle are very similar at the two temperatures. 

Discussion 

Though the precise causes are not fully understood, it is generally agreed that the environmental 
factors of a given experimental setup affect the conductance through a molecule. The conduc- 
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tance of a bridged molecule diluted in a monolayer of non-bridging adsorbate molecules has been 
shown to change when different adsorbate molecules are employed;'^ this result was explained 
by changes in the relative surface coverage for different adsorbates, which can alter the electrode 
work function. Building from this body of work, our results suggest that changes in the electrode 
work function may not be the only factor affecting conductance, as the bonding geometry and tilt 
angle of bridged molecules are both influenced by monolayer density. In particular, monolayer 
density influences whether a molecule is able to sample the specific geometries required for bridg- 
ing while also affecting the availability of preferred adsorption sites. The detailed atomic structure 
of the electrodes also influences the availability of bonding sites. Note that electrode geometry 
and bonding geometry have been linked previously .122 Haiss et al.^^ performed measurements of 
single-molecule conductance using four different experimental techniques, with each method pro- 
ducing different relative populations of the conductance histogram peaks. The authors ascribed this 
behavior to changes in the electrode geometry between methods, which affected the most probable 
bonding geometries. Our results also indicate that electrode geometry affects the most probable 
bonding geometries. While ideal, atomically sharp tips predict a predominance of on-top/on-top 
bonding geometry at high Z, ruptured NW tip geometries allow on-bridge bonding at high Z. 

However, environmental factors do not greatly affect the properties of a junction in all cases. 
For instance, with BDT it may be reasonable to ignore certain environmental effects for inter- 
electrode separations of greater than ~9.5 A. In this regime, the tilt angle is similar regardless of the 
temperature and whether adsorbate molecules are present, and the probability of forming a multi- 
molecule junction is low; this makes Z > 9.5 A well-suited for comparisons between experiment, 
theory, and simulation, since simplified treatments of the junction environment do not significantly 
affect the properties. On the other hand, for Z < 9.5 A, tilt angle data diverge in cases where 
monolayer effects are ignored and the probability of forming multi-molecule junctions increases 
appreciably. In this regime, using simplified treatments for the junction environment may result in 
inaccurate predictions of structure, and thus give rise to incorrect conductance results. In this case 
it is necessary to perform simulations, such as those we describe in this paper, to provide input or 
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guidance for determining the most probable structures for theoretical calculations. 

Although our simulations more closely resemble the MCBJ experimental technique than pre- 
vious simulation studies, we note a few important differences between our method and the exper- 
iments. The first difference is that unlike MCBJ experiments, the electrodes in our simulations 
are not contacted prior to forming each molecular junction. While our simulations often result in 
molecular junctions immediately following NW rupture, for the purposes of gathering meaningful 
statistics and reducing computational expense, we have chosen to simulate the spontaneous forma- 
tion of molecular junctions without contacting the electrodes. In this respect, the junction forma- 
tion process in our simulations is more similar to that of the I{s) and I(t) experimental methods of 
Haiss and co-workers.^^^^ We note that contacting the electrodes may help overcome activation 
barriers involved in junction formation, especially at lower temperature, where the spontaneous 
formation of a molecular junction is less likely (see Figure 5a). 

Another important difference is that while we simulate the compression of a junction, in MCBJ 
experiments the conductance is typically monitored as a junction is elongated. This is an impor- 
tant difference considering the strength of the S-Au bond is high enough to pull short monatomic 
chains of Au atoms out of a surface during elongation,'^2llQ] ^j^^^ ^j^yg j^^y j-ggy^j different elec- 
trode structures than those we use. Despite not considering such effects, we argue that the com- 
pression of a junction prior to electrode contact is a fundamental aspect of the experiments that 
is likely to influence the structures emerging during elongation. Therefore, investigating the de- 
tails of the compression process is essential to understanding the behavior of molecular junctions. 
Furthermore, for the results presented here, we fix the structure of the Au tips during the compres- 
sion/bridging portion of our simulation methodology. This appears to be a reasonable assumption, 
as we did not observe significant rearrangements of the tips during test calculations that allowed 
the tip structure to change during compression/bridging. However, it is important to note that ex- 
periments typically span significantly longer timescales than accessible to simulation, thus atomic 
rearrangements of the Au tips may additionally be important. 

Finally, we address the impact of elongation rate. We perform simulations with a fixed elon- 
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gation rate of 1 m/s and temperature of 298 K. The rate of elongation will influence the resulting 
structural evolution of the wire, however, this effect will be much more significant at low temper- 
ature (e.g., ~4 K), in accordance with the universal energy release mechanism.l^^^ Moreover, in 
our previous simulations of Au NWs elongated at 298 K in vacuum, we did not observe signifi- 
cant differences in the spectrum of resulting tip geometries for rates ranging from 0.033 to 2 m/s, 
although subtle differences in the elongation pathway were observed.*^ We have chosen 1 m/s for 
this study as it is more computationally tractable for systems including BDT. 

Conclusion 

We have demonstrated the utility of a new simulation method that allows for the incorporation 
of important environmental factors into simulations of the formation of molecular junctions. Us- 
ing this tool, we studied aspects of molecular junctions previously inaccessible with simulation. 
We showed that the extent of surface coverage affects the number of bridged molecules. Single- 
molecule junctions were found to occur commonly at intermediate to high surface coverages; 
however, at low inter-electrode separations maximum surface coverage was found to provide the 
highest probability of yielding single-molecule junctions, owing to the limited occurrence of multi- 
molecule junctions in densely packed monolayers. We found that the reduced adsorbate-adsorbate 
interactions at low to intermediate surface coverages leads to relatively high probabilities for form- 
ing multi-molecule junctions. We demonstrated that electrode geometry affects the number, bond- 
ing geometry, and tilt angle of bridged molecules. In addition to influencing the number of bridged 
molecules, monolayer interactions were found to give rise to bonding geometry that is higher in 
energy than the preferred bonding geometry and tilt angles that are higher than those of bridged 
molecules in absence of a monolayer. These are important findings since it has been previously 
demonstrated that both bonding geometry and tilt angle can affect conductance by at least an order 
of magnitude,! ^ ^ I ^'^ l '^^l ^ while also impacting the measured lETS."^ Finally, we demonstrated that 
a low temperature (77 K) significantly reduces the number of bridged molecules, while resulting 
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in only small changes in the bonding geometry and tilt angle, in comparison to 298 K. Our results 
offer guidance on the design of monolayers and electrode geometries to yield desired properties, 
such as specific bonding geometries and/or tilt angles to control conductance. While in this paper 
the focus was on Au-BDT systems, we emphasize that the simulation methodology is in principle 
applicable to any metal-molecule system. 

Methods 

Elongation and Rupture of BDT- Coated Au NWs 

In the hybrid MD-MC scheme, two processes are simulated that are not part of a typical MD sim- 
ulation: elongation of a Au NW and chemisorption of BDT molecules onto a Au NW surface. We 
implement Au NW elongation in the presence of BDT by performing MD simulations within the 
LAMMPS simulation package in the canonical ensemble (constant NVT) with the Nose-Hoover 
thermostat applied to control the temperature at 298 K. The equations of motion are integrated 
using the velocity Verlet algorithm and rRESPA multi-timescale integrator, with outer and inner 
loop time steps of 2.0 and 0.4 fs, respectively. Elongation is carried out using a stretch-and-relax 
technique in which two layers of rigid "gripping" atoms on one end of the wire are periodically 
displaced a small amount (0.1 A) in the [001] direction. Two additional layers of rigid "grip- 
ping" atoms reside on the opposite end of the wire, while atoms within the core of the wire are 
dynamic. We allow the core atoms to relax for 10 ps between displacements of the "gripping" 
atoms, corresponding to a nominal elongation rate of 1.0 m/s. Au-Au interactions are modeled us- 
ing the second-moment approximation to the tight-binding potential (TB-SMA).'^TB-SMA is a 
semi-empirical many-body potential capable of capturing the band character of metallic elements 
at a relatively low computational cost. Furthermore, when compared to DFT calculations, TB- 
SMA outperforms other many-body potentials at describing the energetic and structural evolution 
of elongating Au NWs.l^BDT bonding to the Au NW is held fixed during MD runs; that is, the 
site at which a BDT molecule is bonded to the Au NW does not change. 
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In order to model BDT chemisorption, every ten elongations (i.e., every 1 A of NW elongation) 
we perform MC sampling, with 60,000 fixed-jLiVr (where ji is the chemical potential of a BDT 
molecule, V is volume, and T is temperature) moves followed by 160,000 fixed-A^Vr (where 
is the number of BDT molecules) moves. We find that applying this MC protocol is sufficient for 
equilibration of the metal- molecule interface. That is, applying more MC moves and/or applying 
MC moves at more frequent elongation intervals does not change the results significantly. The 
computational cost is also reasonable with this MC sampling frequency. 

The purpose of the MC simulations is to allow the BDT to efficiently sample favorable bond- 
ing sites on the Au NW surface. Additionally, since the number of BDT molecules is allowed 
to change during fixed- jlVT MC moves, the density of BDT in the simulation box remains rela- 
tively constant for the duration of the elongation/rupture process. During fixed-/iVr MC moves, 
we select a given move type with probabilities 0.45, 0.45, 0.04, 0.04, and 0.02 for BDT center- 
of-mass (COM) displacement, COM rotation, insertion, deletion, and identity swap, respectively. 
For fixed-NVT MC moves we select move types with probabilities 0.49, 0.49, and 0.02 for BDT 
COM displacement, COM rotation, and identity swap, respectively. In accordance with the pre- 
vious work of Pu et a/.'^the excess chemical potential, jigx, of both BDT species is set to -0.525 
kcal/mol. 

We generate the initial wire configuration by taking a cylindrical cut from a FCC lattice, ori- 
ented along the [001] direction. The NW contains a total of 2070 Au atoms, and is long enough 
(12.3 nm) to avoid boundary effects. The box dimensions in the x and y directions are set to 5 nm, 
while the length of the box in the z direction is variable. We apply periodic boundary conditions in 
all three directions. 

Formation of Molecular Junctions 

Following NW rupture, each BDT-functionalized gold tip is allowed to relax its structure at 298 K 
using MD. Since molecular junctions form locally in the break junction created by NW rupture, 
we extract 100 Au atoms at each tip prior to pushing the tips together, which considerably reduces 
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the computational rigor of the simulations. The target surface coverage is obtained by performing 
MC simulations at constant jiVT. To obtain different monolayer arrangements on the tips, we 
remove all chemisorbed BDT and perform the simulations with the bare Au tip as the starting 
point, initializing the psuedorandom number generator of each simulation with a different random 
seed. Next, the bulk BDT is "evaporated" from the simulation box, which is a standard practice in 
real experiments. In the simulation, "evaporation" is accomplished simply by removing from the 
simulation box all of the BDT molecules not bonded to one or more Au atoms. Following this, 
the BDT SAM is equilibrated at constant NVT for 20 million MC moves. Only one S atom in a 
BDT molecule is allowed to bond to the electrode during this process; however, BDT molecules 
are allowed to lie flat on an electrode with both S atoms bonded during the subsequent molecular 
junction formation runs. In all cases the maximum BDT displacement and rotation is adjusted to 
obtain a 40% acceptance rate. The Au atoms are held fixed while the BDT molecules are modeled 
as rigid molecules from an optimized structure using the universal force field.!^ Periodic boundary 
conditions are applied in the x and y directions, while reflective walls are placed parallel to the x-y 
plane at the +z and -z boundaries. We use a box size of 3.5 nm in the x and y directions. 

The spontaneous formation of a molecular junction at fixed inter-electrode distance occurs on 
time scales of ~0.1 s in experiment; this includes time required for bond formation and for 
the molecule to explore sufficient phase space for bridging. Unfortunately, these time scales are 
inaccessible with molecular dynamics simulations, where time steps for integrating the equations 
of motion are typically on the order of 10^^^ s. To access time scales of experiments, we model 
the formation of a S-Au bond using a bonding cutoff, such that if a S atom moves within 3.66 A of 
the appropriate bonding site (on-top or on-bridge, depending on the identity of the molecule), the 
S-H bond dissociates (with the H atom discarded from the simulation) and the S atom covalently 
bonds to the Au site. The MC method does not provide information about the dynamics of bond 
formation, but rather produces thermodynamically favored, equilibrium configurations. We expect 
this approximation to be valid in the limit of slow compression rates, where experimental systems 
are given sufficient time to reach equilibrium. We note that this treatment may slightly overpredict 
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the formation of molecular junctions since we do not consider the details of the BDT chemisorption 
process, which is beyond the scope of conventional CFF methods; however, we expect the trends 
to be qualitatively valid, as our treatment of the bonding process is consistent for all systems. 
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